Image processing method and shading reference data generation method

ABSTRACT

A technology capable of effectively and stably correcting brightness unevenness in an image obtained by imaging is provided. An image processing method includes an image obtaining step of obtaining an original image obtained by imaging an imaging object together with a substantially uniform background, an approximation step of specifying an approximation formula for approximating a two-dimensional luminance profile of the original image to a probability density function of a two-dimensional normal distribution, and a correction step of correcting a luminance value of each pixel constituting the original image on the basis of a luminance value of the two-dimensional luminance profile expressed by the approximation formula.

TECHNICAL FIELD

This invention relates to a technology for solving brightness unevennessby correcting an image obtained by imaging.

BACKGROUND

In a technology for imaging an imaging object under illumination,brightness unevenness may appear in an image due to properties of anillumination system or an imaging system. Even if a background of animaging object has a substantially uniform optical density, suchbrightness unevenness possibly occurs, for example, due to lightquantity unevenness or the like of illumination light incident on theimaging object. A shading correction technology is known as a technologyfor correcting such brightness unevenness. For example, in a technologydescribed in patent literature 1, shading correction data generated onthe basis of shading properties of an image imaged without a specimen isapplied to each pixel of an image obtained by imaging the specimen(imaging object) using a microscope. In this way, image densityunevenness due to an optical system is solved.

CITATION LIST Patent Literature

[Patent literature 1] JP 2005-072967A

SUMMARY Technical Problem

For example, in the case of imaging an imaging object in a liquid,illumination light is incident on the imaging object via a liquidsurface, but a state of the liquid surface is not necessarily constant.As in this example, incident conditions of illumination light may bedifferent for each imaging object. In such a case, the shadingproperties obtained beforehand are not necessarily suitable and aneffective correction may not be possible.

Concerning this problem, in the conventional technology described above,an operator can adjust correction properties while confirming an imageafter correction displayed on a monitor. However, expert knowledge isnecessary to efficiently perform an adjusting operation and a correctionresult depends on the operator's subjectivity. Thus, there is a problemin the stability of the correction result.

Solution to Problem

This invention was developed in view of the above problem and aims toprovide a technology capable of effectively and stably correctingbrightness unevenness in an image obtained by imaging.

To achieve the above object, an aspect of the present invention is animage processing method which includes: an image obtaining step ofobtaining an original image obtained by imaging an imaging objecttogether with a substantially uniform background; an approximation stepof specifying an approximation formula for approximating atwo-dimensional luminance profile of the original image to a probabilitydensity function of a two-dimensional normal distribution; and acorrection step of correcting a luminance value of each pixelconstituting the original image on the basis of a luminance value of thetwo-dimensional luminance profile expressed by the approximationformula.

Another aspect of the present invention is a shading reference datageneration method for generating shading reference data serving as areference when a shading correction is performed on an original imageobtained by imaging an imaging object together with a substantiallyuniform background. To achieve the above object, the shading referencedata generation method includes: an image obtaining step of obtainingthe original image; an approximation step of specifying an approximationformula for approximating a two-dimensional luminance profile of theoriginal image to a probability density function of a two-dimensionalnormal distribution; and a setting step of setting data corresponding toa virtual image having the two-dimensional luminance profile expressedby the approximation formula as the shading reference data.

In the invention thus configured, the approximation formula of theprobability density function expressing the two-dimensional normaldistribution is specified on the basis of the two-dimensional luminanceprofile of the original image under the assumption that a luminanceprofile of a background part of the imaging object follows thetwo-dimensional normal distribution. The two-dimensional luminanceprofile expressed by this approximation formula has a meaning as amathematical model approximately expressing shading properties.According to the knowledge of the inventors of this application, suchapproximation is effective in many imaging situations where a backgroundis substantially uniform, for example, in the imaging of an imagingobject in a liquid. Specifically, a correction similar to that of theshading correction technology is performed using the luminance profileexpressed by the approximation formula obtained as described above inthe invention instead of a profile corresponding to shading propertiesof an optical system obtained by imaging beforehand in the conventionaltechnology. In this way, brightness unevenness in the original image canbe solved and an image having uniformized brightness can be obtained.Since a process does not depend on an operator's subjective judgment, astable result can be obtained.

Advantageous Effects of Invention

As described above, according to the image processing method of theinvention, brightness unevenness in an image obtained by imaging can beeffectively and stably corrected. Further, according to the shadingreference data generation method of the invention, shading referencedata suitable for correcting brightness unevenness in an image can begenerated.

The above and further objects and novel features of the invention willmore fully appear from the following detailed description when the sameis read in connection with the accompanying drawing. It is to beexpressly understood, however, that the drawing is for purpose ofillustration only and is not intended as a definition of the limits ofthe invention.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a diagram showing a schematic configuration of an imagingapparatus to which an imaging method according to the invention can beapplied.

FIG. 2 is a diagram showing an example of a luminance profile of a wellimage.

FIGS. 3A and 3B are graphs showing the probability density function ofthe two-dimensional normal distribution.

FIG. 4 is a flow chart showing the shading correction process in thisembodiment.

FIG. 5 is a flow chart showing an example of the modeling process.

FIGS. 6A to 6C are diagrams showing a method for obtaining initialvalues of the recurrence formula from the process target image.

FIGS. 7A and 7B are diagrams showing an example of a binarized image.

FIGS. 8A and 8B are diagrams showing an example of a corrected image.

DESCRIPTION OF EMBODIMENTS

FIG. 1 is a diagram showing a schematic configuration of an imagingapparatus to which an imaging method according to the invention can beapplied. This imaging apparatus 1 is an apparatus for imaging a cell,cell colony, microbe or the like (hereinafter, called generically “cellor the like”) cultured in liquid poured into recesses called wells Wformed on the upper surface of a well plate WP.

The well plate WP is generally used in the fields of drug discovery andbioscience. A plurality of wells W having a substantially circularcross-section and a transparent and flat bottom surface are disposed tothe upper surface of a plate having a flat plate shape. The number ofthe wells W on the well plate WP is arbitrary. For example, a well plateWP having 96 (12×8 matrix array) wells can be used. A diameter and adepth of each well W are typically about several mm. Note that the sizeof a well plate and the number of wells used in this imaging apparatus 1are arbitrary without being limited to these. For example, well platehaving 6 through 384 wells may be used. Further, the imaging apparatus 1can be applied to image the cell or the like not only in the well plateprovided with a multitude of wells but also, in a flat container calleda “dish”, for example.

A predetermined amount of liquid as a culture medium M is poured intoeach well of the well plate WP. Cells or the like cultured underpredetermined culture conditions in this liquid are imaging objects ofthis imaging apparatus 1. The culture medium M may be added withappropriate reagents or may be gelled after being poured into the wellsW in a liquid state. In this imaging apparatus 1, for example, a cell orthe like cultured on an inner bottom surface of the well can be animaging object as described later. About 50 to 200 microliters of theliquid is generally usually used.

The imaging apparatus 1 includes a holder 11 which holds the well plateWP, an illuminator 12 arranged above the holder 11, an imager 13arranged below the holder 11 and a controller 14 which includes a CPU141 controlling the operation of these components. The holder 11 holdsthe well plate WP in a substantially horizontal posture by being held incontact with a peripheral edge part of the lower surface of the wellplate WP carrying sample together with liquid in each well W.

The illuminator 12 emits an illumination light toward the well plate WPheld by the holder 11. For example, a white LED (light emitting diode)may be used as a light source of the illumination light. A combinationof the light source and an appropriate illumination optical system areused as the illuminator 12. The cells or the like in the well W disposedto the well plate WP are illuminated by the illuminator 12 from above.

The imager 13 is provided below the well plate WP held by the holder 11.In the imager 13, an imaging optical system is arranged at a positionright below the well plate WP. An optical axis of the imaging opticalsystem extends in a vertical direction. FIG. 1 shows a side view. Aright and left direction of the figure indicates a horizontal directionand an up and down direction of the figure indicates a verticaldirection.

The imaging of the cell or the like in the well W is performed by theimager 13. Specifically, light emitted from the illuminator 12 andincident on the surface of the liquid from above the well W illuminatesthe imaging object. Light transmitted downward from the bottom surfaceof the well W is incident to a light receiving surface of an imagingelement 132 via the imaging optical system including an objective lens131. An image of the imaging object is formed on the light receivingsurface of the imaging element 132 by the imaging optical system isimaged by the imaging element 132. The imaging element 132 is a linearimage sensor having a one-dimensional light receiving surface or an areaimage sensor having a two-dimensional light receiving surface. A CCDsensor or a CMOS sensor can be used as the imaging element 132.

The imager 13 is capable of moving in the horizontal direction and thevertical direction by a mechanism controller 146 provided in thecontroller 14. Specifically, the mechanism controller 146 moves theimager 13 in the horizontal direction by operating a drive mechanism 15based on a control command from the CPU 141. By doing so, the imager 13moves relative to the well W in the horizontal direction. Further,focusing is performed by moving the imager 13 in the vertical direction.

Further, the as indicated by arrows with dotted lines shown in FIG. 1,the driving mechanism 15 moves the illuminator 12 integrally with theimager 13 when the imager 13 is moved in the horizontal direction.Specifically, the illuminator 12 is arranged such that a center ofemitted light substantially coincides with the optical axis of theimaging optical system. When the imager 13 moves in the horizontaldirection, the illuminator 12 also moves in conjunction with the imager13.

When an image of one well W in detail is required, the imaging isperformed in a state that a whole of the well W is included in a fieldof view. In this time, the mechanism controller 146 positions the imager13 in the horizontal direction such that the optical axis of the imagingoptical system coincides with the center of the well W. By doing so,whichever well W is imaged, the center of the well W and the center ofemitted light are always position on the optical axis of the imagingoptical system. Consequently, the illuminating condition becomesconstant regardless of which well W is to be imaged, wherefore imagingconditions can be maintained to be satisfactory.

On the other hand, when imaging is required to perform with a pluralityof the wells W more quickly, the imaging is performed with moving theimager 13 horizontally. Specifically, the mechanism controller 146realizes a scan movement by moving the illuminator 12 and the imager 13integrally with constant speed in the horizontal direction. The imager13 performs imaging periodically, thereby a wider area is imaged.Particularly, when the imaging device of the imager 13 is theone-dimensional image sensor, a two-dimensional image can be obtained bythe scan movement of the imager 13 relative to the imaging object.

The image signal output from the imaging device of the imager 13 is sendto the controller 14. Specifically, the image signal is input to an ADconverter (A/D) 143 provided in the controller 14 and converted intodigital image data. The CPU 141 performs appropriate image processingsbased on the received image data. The controller 14 further includes animage memory 144 for storing image data and a memory 145 for storingprograms to be executed by the CPU 141 and data generated by the CPU141, but these may be integrated. The CPU 141 performs variablecalculation processings described later by executing a control programstored in the memory 145.

Besides, the controller 14 is provided with an interface (I/F) 142. Theinterface 142 has a function of performing data exchange with anexternal apparatus connected via a communication line besides a functionof receiving an operation input from a user and presenting informationsuch as processing results to the user. To realize the user interfacefunction, an input receiver 147 for receiving an operation input fromthe user and a display 148 for displaying the messages to the user, aprocessing result or the like are connected to the interface 142.

Note that the controller 14 may be an exclusive device including abovehardware or may be a general-purpose processing device such as apersonal computer or a workstation installed with the control programfor performing the process described above. Specifically, ageneral-purpose computer apparatus may be used as the controller 14 ofthe imaging apparatus 1. When a general-purpose processing device isused as the controller 14, the imaging apparatus 1 may have just aminimal control function for controlling each components of the imager13.

FIG. 2 is a diagram showing an example of a luminance profile of a wellimage. As shown in an upper part of FIG. 2, a culture medium M isinjected into a well W to be imaged, cells or the like C, which areimaging objects, are distributed inside the culture medium M, generallyon an inner bottom surface Wb of the well W. The culture medium M is aliquid or a gel obtained by gelling an injected liquid, and the liquidsurface thereof is not necessarily a horizontal surface and generallyforms a meniscus concave upward.

Thus, illumination light Li incident on the liquid surface from abovesubstantially propagates straight in a part near a center of the well Wand transmits to a side below the well W via the well inner bottomsurface Wb. In contrast, the illumination light Li is refracted at theliquid surface at a position near a side wall surface of the well W anda propagation direction thereof is changed to an outward direction fromthe center of the well W. As a result, in an imaged image, a relativelyhigh luminance is obtained near the center of the well W, whereas theluminance decreases toward a peripheral edge part of the well W as shownin an example of a luminance profile in a lower part of FIG. 2. Besides,wiggling variations of the luminance due to the cells or the like C inthe culture medium M, noise or the like possibly appear in the image.

The culture medium M itself is substantially homogenous and an opticaldensity thereof is thought to be roughly uniform. Thus, an image densityof the culture medium M serving as a background of the cells or the likeC is also supposed to be inherently uniform. However, due to theaforementioned reason, brightness unevenness of, for example, beingbright in a central part and dark in a peripheral edge part occurs in anactual image. Further, brightness unevenness possibly occurs also due tothe inclination of an incident direction of illumination light.

A technology for correcting such brightness unevenness is a shadingcorrection technology. However, the amount and surface state of theculture medium M are not necessarily constant such as due to variationswhen the liquid is injected, and an incident light quantity distributionof illumination light may differ for each well W. Due to these causes,it may not be suitable in some cases to apply a general shadingcorrection technology. Specifically, a reference image obtained underillumination conditions serving as a reference is necessary in aconventional shading correction technology, but the reference cannot beuniquely determined in the above cases.

Accordingly, in this embodiment, data corresponding to a reference imageis extracted to perform a correction from an original image, using animage obtained by imaging the well W containing the imaging objects asthe original image. Specifically, assuming that brightness unevennessdue to the surface state of the culture medium M and the illuminationlight quantity distribution has a luminance distribution according to atwo-dimensional normal distribution, a probability density functionapproximately expressing this luminance distribution is specified. Then,the original image is corrected, using a virtual model image having aluminance distribution expressed by the obtained probability densityfunction instead of the reference image. By doing so, an image havingbrightness unevenness in the original image solved is obtained.

FIGS. 3A and 3B are graphs showing the probability density function ofthe two-dimensional normal distribution. More specifically, FIG. 3A is agraph showing a luminance profile according to the two-dimensionalnormal distribution as a three-dimensional map. FIG. 3B is a graphshowing the luminance profile as a two-dimensional map by equi-luminancelines. It is assumed that a coordinate of the image in a horizontaldirection is an X coordinate and a coordinate thereof in a verticaldirection is a Y coordinate. At this time, a luminance value L(x, y) ofa pixel P(x, y) having an X coordinate value of x and a Y coordinatevalue of y, out of each pixel constituting the image, can be expressedby the following (Equation 1), which is a probability density functioncorresponding to the two-dimensional normal distribution.

F(x,y;a,b,c,d,e,f,g)=a exp{−b(x−e)² +c(y−f)² +d(x−e)(y−f)}+g  (Equation1)

The luminance value L(x, y) can take various values depending on valuesof parameters a, b, c, d, e, f and g in the above (Equation 1), but theluminance profile on an arbitrary straight line in an XY plane is anormal distribution curve as is clear from graph shapes of FIGS. 3A and3B. Further, an equi-luminance plane generally has an elliptical shape,and a major axis direction thereof can be any of various directions.

If (Equation 1) and a curved surface shown in FIG. 3A are compared, theparameter a corresponds to a luminance value at a peak of the curvedsurface, and the parameters b, c respectively correspond to spreadshapes of the curved surface in X and Y directions. Further, theparameters e, f correspond to coordinate positions of the peak, and theparameter d represents the rotation of an ellipse main axis. Further,the parameter g corresponds to a luminance value to which the base ofthe curved surface comes closer.

The least-squares method is widely used as a method for obtaining amathematical expression approximately expressing a distribution ofactual data. However, if approximate representation by a general linearhigh-order polynomial is allowed only for the purpose of simplyminimizing an error, it may be possibly not suitable to express aluminance distribution in a well in some cases. More specifically, in abackground part of the image obtained by imaging the inside of the wellW, the luminance value is highest near the center of the well W and theluminance value converges to a fixed value in a peripheral end part.Unless the luminance values in the central part and the peripheral endpart are reflected as constraint conditions in an approximation formula,an image after correction may become unnatural.

In this point, it is reasonable to use the probability density functionof the two-dimensional normal distribution expressed by (Equation 1) asthe approximation formula. Further, the inventors of this applicationevaluated a correction method in this embodiment using images obtainedby imaging various specimens and confirmed that a good result wasobtained in many cases.

By approximating brightness unevenness appearing in the original imagedue to imaging conditions to such a luminance profile and performing acorrection process to remove a contribution of this luminance profilefrom the original image, it is possible to obtain an image having anoptical density inherently possessed by the imaging object reflectedtherein. Such a correction process corresponds to a shading correctionprocess based on the reference image obtained beforehand in theconventional technology. More specifically, this correction process is acorrection process using a model image generated on the basis of theoriginal image instead of the reference image used in the shadingcorrection process of the conventional technology. The correctionprocess of this embodiment is referred to as a “shading correctionprocess” in this embodiment for the sake of convenience below.

FIG. 4 is a flow chart showing the shading correction process in thisembodiment. This process is realized by the CPU 141 of the controller 14executing the control program stored in advance in the memory 145.First, an original image is obtained by imaging the well W provided inthe well plate WP together with the contents thereof (Step S101). Animage of the well W imaged by the imager 13 can be used as the originalimage. Besides, it is also possible to use an image given from anexternal imaging device or data storage device via the interface unit142. If the obtained image includes a plurality of wells W or includesmany margins other than the well W, a rectangular region including oneentire well W and a minimum margin of this image is cut out as a wellregion and the cut-out partial image is set as an original image (StepS102).

If the original image obtained in this way includes unnecessary objectsclearly different from the imaging object such as air bubbles andforeign matters in the culture medium M, these are removed by anappropriate image processing. In this way, a process target imageserving as a target of the shading correction process in this embodimentis generated (Step S103). In the image processing in this case, knownvarious image processing algorithms called “outlier removing” can be,for example, applied. Further, a process of removing an objectdesignated by a user may be applied. By performing such a process, imagequality after the shading correction process can be improved, but thereare many cases where an image includes no such unnecessary object. Thus,this process may be omitted.

Further, the original image also includes the cells or the like servingas the imaging objects. Out of those, objects taking up a particularlylarge area in the well may be removed similarly to the unnecessaryobjects. However, if the imaging objects are general cells, a densitydifference from the culture medium M is small. Thus, even if the imagingobjects remain as they are, an influence on a processing result isthought to be small. This is because a small density change due to thecells or the like is substantially ignored since the process isperformed, regarding the luminance distribution of the background partin the well as the two-dimensional normal distribution in thisembodiment as described below.

Subsequently, a modeling process for generating the model imageequivalent to the reference image in the conventional shading correctionprocess is performed on the basis of the generated process target image(Step S104). The principle of the modeling process is described below.The modeling process is a process for specifying a mathematical modelrepresenting a luminance distribution of a background part (hereinafter,referred to as a “background model”), regarding that the well regioncorresponding to the well W in the process target image is a regionwhere the cells or the like as the imaging objects are distributed inthe background having a predetermined luminance distribution. Asdescribed above, in this embodiment, it is assumed that the luminancedistribution of the background part follows the two-dimensional normaldistribution expressed by (Equation 1). Thus, to specify the backgroundmodel, the unknown parameters a, b, c, d, e, f and g in (Equation 1) maybe determined to fit well to a luminance distribution in the processtarget image.

Fitting is performed by determining the unknown parameters a, b, c, d,e, f and g to minimize a square error of a luminance value L(x, y) ofeach pixel P(x, y) of the process target image and a luminance valueF(x, y) of a corresponding pixel in the background model. The squareerror E can be expressed by the following (Equation 2).

$\begin{matrix}{\mspace{79mu} {{E = {\sum\limits_{i = 1}^{n}\; \left\{ {{F\left( {x_{i},y_{i},{;a_{0}},b_{0},c_{0},d_{0},e_{0},f_{0},g_{0}} \right)} - {\Delta \; F}} \right\}^{2}}}{{{where}\mspace{14mu} \Delta \; F} = {{\frac{\partial F}{\partial a}\Delta \; a} + {\frac{\partial F}{\partial b}\Delta \; b} + {\frac{\partial F}{\partial c}\Delta \; c} + {\frac{\partial F}{\partial d}\Delta \; d} + {\frac{\partial F}{\partial e}\Delta \; e} + {\frac{\partial F}{\partial f}\Delta \; f} + {\frac{\partial F}{\partial g}\Delta \; g}}}}} & \left( {{Equation}\mspace{14mu} 2} \right)\end{matrix}$

In (Equation 2), a sign i is a suffix for distinguishing the individualpixels constituting the image. The process target image is composed of npixels and each pixel is distinguished by the suffix i (i=1, 2, . . . ,n). Coordinate values of the i^(th) pixel Pi is expressed by (x_(i),y_(i)). Further, signs a₀ to g₀ denote specific numerical values asinitial values of the respective parameters a to g. Further, signs Δa toΔg represent differences of the respective parameters a to g.

Further, partial differential coefficients included in each term of anequation expressed by a sign ΔF in (Equation 2) can be respectivelyexpressed as follows from (Equation 1).

$\left. \mspace{785mu} {\left( {{Equation}\mspace{14mu} 3} \right)\begin{matrix}{\frac{\partial F}{\partial a} = {\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}} \\{\frac{\partial F}{\partial b} = {{- {a\left( {x - e} \right)}^{2}}{\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}}} \\{\frac{\partial F}{\partial c} = {{- {a\left( {y - f} \right)}^{2}}{\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}}} \\{\frac{\partial F}{\partial d} = {{a\left( {x - e} \right)}\left( {y - f} \right){\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}}} \\{\frac{\partial F}{\partial e} = {a\left\{ {{2\; {{be}\left( {x - e} \right)}} - {d\left( {y - f} \right)}} \right\}^{2}{\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}}} \\{\frac{\partial F}{\partial f} = {a\left\{ {{2\; {{cf}\left( {y - f} \right)}} - {d\left( {x - e} \right)}} \right\}^{2}{\exp \left( {{- {b\left( {x - e} \right)}^{2}} - {c\left( {y - f} \right)}^{2} + {{d\left( {x - e} \right)}\left( {y - f} \right)}} \right)}}}\end{matrix}} \right\}$

The following equation is obtained from a condition that the right sideof (Equation 2) becomes zero by being partially differentiated by Δa toΔg.

                                             (Equation  4)${\begin{pmatrix}\left\lbrack \left( \frac{\partial F}{\partial a} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial a}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial b} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial b}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial c} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial c}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial d} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial d}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial e} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial e}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial f} \right)^{2} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial f}\frac{\partial F}{\partial g}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial a}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial b}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial c}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial d}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial e}} \right\rbrack & \left\lbrack {\frac{\partial F}{\partial g}\frac{\partial F}{\partial f}} \right\rbrack & \left\lbrack \left( \frac{\partial F}{\partial g} \right)^{2} \right\rbrack\end{pmatrix}\begin{pmatrix}{\Delta \; a_{0}} \\{\Delta \; b_{0}} \\{\Delta \; c_{0}} \\{\Delta \; d_{0}} \\{\Delta \; e_{0}} \\{\Delta \; f_{0}} \\{\Delta \; g_{0}}\end{pmatrix}} = \begin{pmatrix}\left\lbrack {\frac{\partial F}{\partial a}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial b}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial c}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial d}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial e}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial f}F_{i}} \right\rbrack \\\left\lbrack {\frac{\partial F}{\partial g}F_{i}} \right\rbrack\end{pmatrix}$

An operant [Q] written in square brackets in (Equation 4) indicates atotal value of values of the function Q obtained for all the pixels.Further, signs Δa₀ to Δg₀ indicate difference values of the parameters ato g from the respective initial values. Further, on the right side,

F _(i) =F(x _(i) ,y _(i)).

Values of partial differential coefficients in (Equation 4) can beobtained as numerical values by substituting the initial values a₀ to g₀and the coordinate values (x_(i), y_(i)) of each point into (Equation3). Note that, in each component of a 7×7 matrix constituting the leftside of (Equation 4), a multiplication order of primary differentialvalues is interchangeable. Thus, this matrix is a symmetric matrix.

The background model is determined by solving (Equation 4) and obtainingthe parameters a, b, c, d, e, f and g using the initial values a₀ to g₀and the difference values Δa₀ to Δg₀. In this embodiment, theLevenberg-Marquardt method (hereinafter, abbreviated as the “LM method”)is used as a solution of the equation to well match the luminancedistribution of the background model particularly in the well peripheralend part with an actual measurement (i.e. a background luminancedistribution in the process target image) by reflecting the constraintconditions described above. The LM method is an example of a calculationmethod for converging an unknown parameter to an optimal value byiteratively calculating for a recurrence formula in which an appropriateinitial value is given to the parameter. By setting the initial valuereflecting constraint conditions, a convergence result and an actualmeasurement result can be well matched.

The recurrence formula is expressed by the following (Equation 5) usinga suffix K indicating a term number.

U ^((K+1)) u ^((K))−(H _(u) ^((K)) +hD[H _(u) ^((K))])⁻¹∇_(u) J^((K))  (Equation 5)

In (Equation 5), a symbol Hu represents a 7×7 matrix on the left side of(Equation 4) and a symbol ∇uJ represents a vector on the right side of(Equation 4). Further, a symbol D[Hu] represents a function for taking adiagonal term of the matrix Hu, and a symbol h is a coefficient alsocalled a Marquardt parameter. A specific procedure of the modelingprocess for specifying the background model of the process target imageis described below using (Equation 5).

FIG. 5 is a flow chart showing an example of the modeling process. Thisprocess is realized by the CPU 141 of the controller 14 executing thecontrol program stored in advance in the memory 145. First, whichmathematical model is used to represent the background model isdetermined (Step S201). Generally, various mathematical models areconceivable as the background model, but (Equation 1), i.e. theprobability density function of the two-dimensional normal distributionis used as the background model in this embodiment as described above.In other words, the process shown in FIG. 5 can be utilized also whenthe background model is expressed by another mathematical formula.

Then, each partial differential coefficient in the equation expressingthe square error of the process target image and the background model isobtained using the equation of the background model (Step S202). In thisembodiment, each partial differential coefficient is expressed by(Equation 3). These partial differential coefficients are substitutedinto (Equation 5).

Subsequently, the initial values of the parameters to be given to therecurrence formula of (Equation 5) are set (Step S203). The initialvalues of the parameters a, b, c, d, e, f and g are determined asfollows on the basis of the luminance value of each pixel constitutingthe process target image.

FIGS. 6A to 6C are diagrams showing a method for obtaining initialvalues of the recurrence formula from the process target image. Morespecifically, FIG. 6A shows an example of the process target image andFIG. 6B is a histogram showing a luminance value distribution of thepixels constituting this image. Further, FIG. 6C is a table showingvarious pieces of data obtained from the histogram. As shown in FIG. 6A,a case is considered where an image of the well W bright in a centralpart and dark in a peripheral edge part is a process target image Ip.Here, the image Ip obtained by imaging the well W into which the culturemedium M not containing the cells or the like is injected is used. Alarge defective part seen in a left upper part of the well W is causedby air bubbles generated on the liquid surface.

As shown in FIG. 6B, when a distribution of luminance values in therespective pixels constituting a circular part corresponding to the wellW, out of the process target image Ip, is expressed by a histogram, apeak corresponding to a low luminance part of the peripheral edge partof the well W appears at a relatively low luminance and the luminancevalues of the respective pixel in the well W are widely distributed on aside higher than this peak. This indicates that brightness unevenness islarge in the well W. From this histogram, each piece of data shown inFIG. 6C is obtained. Out of these, a value 29 of a mode is set as aninitial value of the parameter g corresponding to the luminance to whichthe base part of the curved surface of FIG. 3A comes closer. Further, avalue 138 obtained by subtracting the mode from a maximum value is setas an initial value of the parameter a corresponding to a peak height ofthe curved surface of FIG. 3A.

The other parameters relate to an elliptical shape represented bycross-sections when the curved surface of FIG. 3A are cut by planesparallel to the XY plane. To obtain initial values of these, the processtarget image Ip is binarized using an appropriate threshold value and anelliptical area is extracted. Here, a half value in the two-dimensionalnormal distribution represented by the curved surface of FIG. 3A, i.e. avalue equivalent to a middle point between the initial values of theparameters a, g set in advance (=(a−g)/2), is used as the thresholdvalue.

FIGS. 7A and 7B are diagrams showing an example of a binarized image.More specifically, FIG. 7A shows an example of a binarized image Ibobtained from the image of FIG. 6A and FIG. 7B is a table showingvarious pieces of data obtained from this binarized image Ib. As shownin FIG. 7A, in the binarized image Ib generated from the process targetimage Ip, a region having a luminance value equal to or larger than thethreshold value, out of the process target image Ip, is extracted. Thus,a high-luminance part in a central part is mainly extracted and theshape of the thus extracted region represents a schematic shape of theelliptical region described above. Thus, the initial values of theremaining parameters can be obtained from the shape of the extractedhigh-luminance region.

Note that there is not necessarily a single high-luminance regionextracted by binarization due to brightness unevenness appearing in theimage. Accordingly, if there are a plurality of extracted high-luminanceregions, attention is paid to the one having a largest area out ofthose. Further, if it is clear in appearance that a plurality ofseparated high-luminance regions form one continuous region as a whole,for example, as in an air bubble part shown in a left upper part of FIG.7A, those high-luminance regions may be regarded as one high-luminanceregion.

A dotted line of FIG. 7A represents a rectangle (circumscribingrectangle) R circumscribing a high-luminance region Rh extracted bybinarization, and FIG. 7B shows vertex coordinate positions and size ofthis circumscribing rectangle R. A width, i.e. a length in the Xdirection, of the circumscribing rectangle R is 450 and this is set asthe initial value of the parameter b. On the other hand, a height, i.e.a length in the Y direction, of the circumscribing rectangle R is 462and this is set as the initial value of the parameter c. Further, an Xcoordinate value of a centroid position of the circumscribing rectangleR is 301, a Y coordinate value thereof is 335 and these are respectivelyset as the initial values of the parameters e, f. The initial value ofthe parameter d relating to the rotation of the elliptical region is setto 0.

The initial values of the respective parameters are organized asfollows.

a: 138 (difference between the maximum value and the mode in thehistogram)

b: 450 (width of the circumscribing rectangle in the binarized image)

c: 462 (height of the circumscribing rectangle in the binarized image)

d: 0 (no rotation is assumed)

e: 301 (X coordinate of the centroid of the circumscribing rectangle inthe binarized image)

f: 335 (Y coordinate of the centroid of the circumscribing rectangle inthe binarized image)

g: 29 (mode)

Using a set (a, b, c, d, e, f, g) of the parameters having these as theinitial values, a coefficient h is set to an appropriate value, e.g.0.0001 (Step S204) and the calculation by the LM method is performed tosolve (Equation 5) (Step S205), whereby the difference values Δa₀, Δb₀,Δc₀, Δd₀, Δe₀, Δf₀ and Δg₀ are respectively obtained.

By adding or subtracting the obtained difference values Δa₀, Δb₀, Δc₀,Δd₀, Δe₀, Δf₀ and Δg₀ to or from the currently set values of theparameters a, b, c, d, e, f and g, a new parameter set (a, b, c, d, e,f, g) is generated (Step S206). If an execution number of the iterativecalculation by the LM method has reached a predetermined value (YES inStep S207), the process is finished.

Unless the execution number of the iterative calculation by the LMmethod has reached the predetermined value (NO in Step S207), adifference between a value u^((K)) of the current recurrence formula anda newly calculated value u^((K+1)) is calculated as a residual error(Step S208). If the residual error is smaller than a value obtained inthe last calculation, the equation is judged to approach convergence(YES in Step S209). In this case, to avoid large variations of theparameters by one calculation, the coefficient h is reduced to 1/10 ofthe current value (Step S210). On the other hand, if the residual erroris larger than the value obtained in the last calculation, the equationis judged to move away from convergence (NO in Step S209) and thecoefficient h is increased to ten-fold of the current value.

By performing the iterative calculation using the LM method in this way,the parameter set (a, b, c, d, e, f, g) is optimized. As a result, thebackground model expressed by (Equation 1) having these parameterssubstituted thereinto is optimized according to the luminancedistribution in the actually obtained process target image. Thebackground model generated by the modeling process described above isregarded as the virtual reference image. By correcting the luminancevalues of the respective pixels constituting the process target imageusing the luminance values of the pixels of the model image, a correctedimage having the shading correction process applied thereto is obtained.

Referring back to FIG. 4, the shading correction process is furtherdescribed. When the modeling process is finished, the luminance value ofeach pixel of the process target image is corrected on the basis of thebackground model (Step S105). More specifically, the model image havingthe luminance profile possessed by the background model expressed by(Equation 1) and the optimized parameter set (a, b, c, d, e, f, g) isvirtually generated. Each pixel constituting the model image has a valueof F(x, y) obtained by substituting the coordinate position (x, y)thereof into optimized (Equation 1) as the luminance value. Byperforming the shading correction process regarding this model image asthe reference image, the corrected image having brightness unevenness inthe process target image solved is generated.

The correction process in this case can be performed on the basis of aprinciple similar to that of the known shading correction process andcan be performed, for example, as follows. Specifically, it is assumedthat L1 denotes a luminance value of an arbitrary target pixel in theprocess target image and L2 denotes a luminance value of a pixel locatedat a coordinate position corresponding to the coordinate position of thetarget pixel in the model image expressed by the optimized backgroundmodel. At this time, a luminance value L3 of the pixel occupying thiscoordinate position in the corrected image is obtained by the followingequation:

L3=A·(L1/L2)+B  (Equation 6),

whereby the luminance value of the pixel after correction is determined.

Symbols A, B in (Equation 6) are actual constants appropriatelydetermined in advance and, out of these, the constant A is not 0. Byarranging each pixel having the luminance value obtained in this way atthe predetermined coordinate position, the corrected image havingbrightness unevenness in the original image solved is obtained.

Note that the process target image serving as a target of correction inthis process has the unnecessary objects in the well removed in StepS103. The image having clearly unnecessarily objects such as air bubblesand foreign matters removed can be directly used as a target of thecorrection process. On the other hand, only for the purpose of simplygenerating a background model, the modeling process can be performed,using an image having objects corresponding to the cells or the like inthe well as the original imaging objects at least partially removed, animage having a filtering process for smoothing an image density changeapplied thereto or the like as the process target image.

In such a case, since part of image information relating to the imagingobject may be lost, it is not appropriate to use the image used in themodeling process as the target of the correction process. Thus, it isdesirable to directly use the original image obtained in Step S101 orS102 as the target of the correction process or to separately prepare animage having only clearly unnecessary objects removed therefrom as theprocess target image and use this image for the correction process.

FIGS. 8A and 8B are diagrams showing an example of a corrected image.More specifically, FIG. 8A shows an example of a corrected image Icobtained by correcting the process target image shown in FIG. 6A by theshading correction process of this embodiment. Further, FIG. 8B is ahistogram showing a distribution of luminance values of pixelsconstituting this image Ic. In comparison to the process target image Ipshown in FIG. 6A, the inside of the well W is entirely brighter in thecorrected image Ic shown in FIG. 8A and a difference in brightnessbetween the central part and the peripheral edge part of the well issmaller. As shown in FIG. 8B, it is understood that a variation of theluminance values of the respective pixels becomes smaller, brightnessunevenness is solved and the brightness of the image is uniformized.

By performing the above process on the image including the cells or thelike as the imaging objects, an image in a state where the cells or thelike are distributed in the well region having the backgroundsubstantially uniformized is obtained as the corrected image. In thisway, it is possible to generate an image having an imaging objectclarified by correcting brightness unevenness even from an originalimage having nonuniform illumination conditions in the well W, forexample, due to the meniscus of the liquid surface.

As described above, in the above embodiment, Steps S101, S103, S104 andS105 of FIG. 4 respectively correspond to an “image obtaining step”, an“object removal step”, “an approximation step” and a “correction step”of the invention.

Further, in the above embodiment, the well plate WP corresponds to a“specimen container”, and each well W corresponds to a “recess” of theinvention. Further, the X coordinate axis and the Y coordinate axis inthe imaged image respectively correspond to a “first coordinate axis”and a “second coordinate axis” of the invention. Further, thehigh-luminance region Rh in the binarized image Ib of FIG. 7Acorresponds to a “main region” of the invention.

Note that the invention is not limited to the above embodiment andvarious changes other than those described above can be made withoutdeparting from the gist of the invention. For example, a function ofcorrecting the original image using the background model generated onthe basis of the information obtained from the original image isprovided in the above embodiment. However, a function of generating thebackground model and a function of performing a correction on the basisof the background model may be separated. For example, the invention canbe applied also to an apparatus having only a function of supplying datanecessary for correction to an external apparatus or the like in chargeof a shading correction function.

In this case, the data relating to the background model, e.g.(Equation 1) or the image data corresponding to the virtual model imagerepresented by the background model corresponds to “shading referencedata” of the invention. Then, this apparatus serves as a main bodyimplementing the “shading reference data generation method” according tothe invention. In such a mode, Step S104 corresponds to a “setting step”of the invention.

Further, in the above embodiment, one entire well W is imaged by oneimaging. On the other hand, in the case of imaging wells (e.g. about 6to 48 wells arranged in one well plate) or dishes having a largerdiameter, a range to be imaged may become larger than a visual field ofthe imager 13. In this case, the range to be imaged is imaged by beingdivided into a plurality of partial images and those partial images arecombined by an image processing to generate an image covering the entirerange to be imaged. Even in such a case, the correction processaccording to the invention can be applied to the individual partialimages. By combining the corrected partial images, it is possible toobtain an image with less brightness unevenness in its entirety.

Further, the above embodiment relates to the apparatus with the imagingfunction of imaging the original image. However, the “image processingmethod” and the “shading reference data generation method” of theinvention can also be realized by an apparatus having no imagingfunction itself and configured to perform a process on the basis oforiginal image data given from an external apparatus or the like.

Further, in the modeling process of the above embodiment, it is assumedthat the background model is expressed by the probability densityfunction of the two-dimensional normal distribution. However, themodeling process itself shown in FIG. 5 can be similarly applied to abackground model expressed by an arbitrary mathematical expression. Forexample, it is possible to prepare a plurality of background modelsincluding the two-dimensional normal distribution model in advance andselect any one of the background models in Step S201.

Further, in the correction process (Step S105 of FIG. 4) of the aboveembodiment, the model image is generated from the background modelexpressed by (Equation 1) with the optimized parameters. However, it issufficient to provide information necessary to correct the luminancevalue to each pixel of the process target image and it is not alwaysnecessary to generate the image data of the virtual model image. Forexample, the coordinate value may be substituted into (Equation 1) foreach target pixel and the luminance value of this pixel in the modelimage may be calculated every time. Further, the value obtained bysubstituting the coordinate value into optimized (Equation 1)corresponds to L2 of (Equation 6). Accordingly, a pixel value aftercorrection may be directly calculated from the luminance value and thecoordinate value of the target pixel, using a calculation formulaobtained, for example, by substituting the right side of (Equation 1)into L2 of (Equation 6).

Further, in the above embodiment, the luminance value after correctionof each pixel is calculated using (Equation 6). However, this method ismerely an example and specific contents of the correction process afterthe background model is specified are not limited to the above, butarbitrary.

Although the invention has been described with reference to specificembodiments, this description is not meant to be construed in a limitingsense. Various modifications of the disclosed embodiment, as well asother embodiments of the present invention, will become apparent topersons skilled in the art upon reference to the description of theinvention. It is therefore contemplated that the appended claims willcover any such modifications or embodiments as fall within the truescope of the invention.

INDUSTRIAL APPLICABILITY

The invention is suitable in technical fields requiring brightnessunevenness in an image obtained by imaging, particularly an imageobtained by imaging an imaging object together with a substantiallyuniform background to be solved.

REFERENCE SIGNS LIST

-   1 imaging apparatus-   12 illuminator-   13 imager-   14 controller-   Rh high-luminance region (main region)-   S101 image obtaining step-   S103 object removal step-   S104 approximation step, setting step-   S105 correction step-   W well (recess)-   WP well plate (specimen container)

1. An image processing method, comprising: an image obtaining step ofobtaining an original image obtained by imaging an imaging objecttogether with a substantially uniform background; an approximation stepof specifying an approximation formula for approximating atwo-dimensional luminance profile of the original image to a probabilitydensity function of a two-dimensional normal distribution; and acorrection step of correcting a luminance value of each pixelconstituting the original image on the basis of a luminance value of thetwo-dimensional luminance profile expressed by the approximationformula.
 2. The image processing method according to claim 1, wherein,in the approximation step, the approximation formula is specified byobtaining unknown parameters a, b, c, d, e, f and g using the followingequation:F(x,y)=a·exp{−b(x−e)² +c(y−f)² +d(x−e)(y−f)}+g as the approximationformula when x denotes a coordinate on a first coordinate axis of theoriginal image and y denotes a coordinate on a second coordinate axisorthogonal to the first coordinate axis.
 3. The image processing methodaccording to claim 2, wherein the approximation formula is specified bythe Levenberg-Marquardt method in the approximation step.
 4. The imageprocessing method according to claim 3, wherein: a maximum value, aminimum value and a mode in a luminance value histogram of the originalimage are obtained in the approximation step; and when the originalimage is binarized using a half value of a difference between themaximum value and the mode as a threshold value and a region having alargest area, out of regions in which pixels having a luminance valueequal to or larger than the threshold value continue, is set as a mainregion, initial values of the respective unknown parameters are set asfollows: a: the difference between the maximum value and the mode, b: alength of the main region in a direction of the first coordinate axis,c: a length of the main region in a direction of the second coordinateaxis, d: 0 e: a coordinate value of a centroid of the main region on thefirst coordinate axis, f: a coordinate value of the centroid of the mainregion on the second coordinate axis, and g: the mode.
 5. The imageprocessing method according to claim 1, wherein, in the correction step,when L1 denotes the luminance value of each pixel of the original imageand L2 denotes a luminance value of a corresponding pixel in thetwo-dimensional luminance profile expressed by the approximationformula, a luminance value L3 of the pixel after correction is obtainedby the following equation:L3=A·(L1/L2)+B using actual constants A (A≠0) and B determined inadvance.
 6. The image processing method according to claim 1, comprisingan object removal step of removing at least one of objects included inthe original image from the original image prior to the approximationstep.
 7. The image processing method according to claim 1, wherein theoriginal image is an image obtained with a biological specimen in aliquid as the imaging object.
 8. The image processing method accordingto claim 7, wherein the original image is an image obtained by imagingthe liquid carried in a recess provided in a specimen container andhaving a substantially circular horizontal cross-sectional shapetogether with the entire recess.
 9. The image processing methodaccording to claim 7, wherein the original image is an image imaged in astate where illumination light is incident on the imaging object via aliquid surface of the liquid.
 10. A shading reference data generationmethod for generating shading reference data serving as a reference whena shading correction is performed on an original image obtained byimaging an imaging object together with a substantially uniformbackground, comprising: an image obtaining step of obtaining theoriginal image; an approximation step of specifying an approximationformula for approximating a two-dimensional luminance profile of theoriginal image to a probability density function of a two-dimensionalnormal distribution; and a setting step of setting data corresponding toa virtual image having the two-dimensional luminance profile expressedby the approximation formula as the shading reference data.
 11. Theshading reference data generation method according to claim 10,comprising an object removal step of removing at least one of objectsincluded in the original image from the original image prior to theapproximation step.